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SYSTEMS 
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' Abstract. We study the convergence of a class of asymptotic preserving numerical schemes initially 

, proposed by F. Filbet & S. Jin [23] and G. Dimarco & L. Pareschi [l7] in the context of nonlinear and stiff 

kinetic equations. Here, our analysis is devoted to the approximation of a system of transport equations 
I with a nonlinear source term, for which the asymptotic limit is given by a conservation laws. We investigate 

the convergence of the approximate solution to a nonlinear relaxation system, where e > is a 

physical parameter and h represents the discretization parameter. Uniform convergence with respect to e 
fT^ , and h is proved and error estimates are also obtained. Finally, several numerical tests are performed to 

illustrate the accuracy and efficiency of such a scheme. 
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O ' 1. Introduction 

Many physical problems are governed by hyperbolic conservation laws with non vanishing stiff source 
terms. These problems can describe the effect of relaxation toward an equilibrium, as in the kinetic 
theory of gases for example. Actually for monatomic gases, when an equilibrium state is perturbed, it 
gradually relaxes to the equilibrium state with Maxwellian velocity distribution. The time scale for such a 
' relaxation process may not be short and the phenomenon of thermo-nonequilibrium becomes important. 
In this case, the compressible Euler or Navier-Stokes equations can be applied. In the domain of kinetic 
equations, the relaxation parameter is represented by the dimensionless Knudsen number e > 0, defined 
as the ratio of the mean free path of the particles over a typical length scale, such as the size of the 
spatial domain. At the continuous level, it has been shown p] that, for the Boltzmann equation 

dtf + v-V^f = ^Q{f), 

when the Knudsen number e goes to zero, the distribution function / converges to a local Maxwellian 
Ai, so the macroscopic model (compressible Euler or Navier-Stokes) becomes more adequate to describe 
the behavior. For more details about fluid dynamic limits of kinetic equation, we refer to [2l[3llll[8]. 

In this context, developing robust numerical schemes for kinetic equations that also work in the fluid 
regime becomes challenging. It has been done in the framework of Asymptotic-Preserving (AP) schemes 
, |27]. As it has already been defined in [27], a numerical scheme is Asymptotic Preserving if 
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(1) it is a suitable scheme for the relaxation problem (the kinetic equation), and when the mesh size 
and time step are fixed, whereas the parameter e goes to zero, the scheme becomes a suitable 
discretization of the limit (equilibrium) problem; 

(2) implicit collision terms can be implemented explicitly. 

Indeed, from a numerical point of view, the treatment of the stiffness can not be done with explicit 
schemes, so mathematicians will rather favor the use of semi-implicit of fully implicit schemes. 

One solution offered by E. Gabetta, L. Pareschi and G. Toscani [25] to design an Asymptotic Preserving 
scheme, was to penalize the nonlinear collision operator Q{f) by a linear function A / and then absorb 
the linearly stiff part into the time variable to remove the stiffness. More recently, F. Filbet & S. Jin |23j 
proposed to penalize the Boltzmann's operator by the BGK operator in order to build stable schemes 
with respect to e > 0. If such schemes are now numerically validated and extensively used to discretize 
kinetic equations [171 [23l [24l [25| 127] . their mathematical study has only been done in some particular 
cases ^26j, because of the complexity of general kinetic equations, where the collision operator is nonlinear. 
In this context, hyperbolic conservation laws represent a simplified framework where a theoretical study 
of relaxation schemes can be done [H [TTl [12l [131 [El [221 [23 [23 [35]. Indeed, there is a strong analogy 
between the local relaxation approximation of conservation laws, initially proposed in [27], and the study 
of fluid dynamical limits of kinetic equations [H]. In the domain of hyperbolic conservation laws with 
stiff source terms, the relaxation parameter plays the role of the Knudsen number in kinetic theory, 
while the source term is the analogous of the collision operator, which we want to be the most general 
possible. Few works are devoted to the mathematical analysis of relaxation scheme for the approximation 
of conservations laws. We refer for instance to the series of paper by D. Aregba-DrioUet and R. Natalini 
[1] and A. Chalabi |1 H 1121 [T3] . But for all of them, the relaxation operator is relatively simple and can 
be easily treated explicitly without any additional computational cost. Moreover, the goal of these works 
mainly consists to approximate the asymptotic limit of the relaxation model and not to develop a scheme 
which is accurate and efficient in various regimes. 

Here we want to focus on the approximation of general relaxation system with a nonlinear and stiff 
source term as in the context of the Boltzmann kinetic equation, which is valid both for rarefied (large 
values of e) and hydrodynamic (e — )• 0) regimes. Therefore, we will consider throughout the rest of the 
paper the following hyperbolic relaxation system. For all (t, x) in x M: 



(1.1) 



r dtu' + d^v" = , 



where a > is a constant coefficient to be discussed later, e is the relaxation parameter and TZ : R x M i— 
M is a nonlinear function. The system is completed with the initial conditions: 



(1.2) 



n^(0,x) = lig(x) , 
u^(0, x) = Vq{x) . 



The system of equations (|l.ip - (|1.2p is often referred as a two velocity kinetic equation. Here we will 
assume that the function TZ £ C^(MxM,M) possesses a unique local equilibrium, restricted to the manifold 
{v = A{u)}, that is, 

(1.3) n{u,v) =0 v = A{u) , 

where j4 is a locally Lipschitz continuous function with ^(0) = 0. Therefore, under some assumptions 
on TZ and on the initial data, the solution (u^,t>^) to (jl.ip - (|1.2p converges to {u,v) with v = A{u) and u 
solution to the conservation laws |141 [M] 



(1.4) 



dtu + dxA{u) = , in 

u{t = 0) = Uq , 
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where the initial datum uq is given by 
(1.5) UQ 

Concerning the mathematical study of the system (ll.ip - (11.2p . we refer for instance to [38] 



limun . 

e-s>0 " 



Theorem 1.1. Assume that the initial datum {uq,Vq) is bounded independently of e in L°° n Sy(M). 
Consider 7^ G C^(M x M,M), which satisfies and VTZ is uniformly hounded with respect to v and 

locally hounded with respect to u such that, for any {u,v) in [—Uq^ Uq] x M; 

f \dun{u,v)\ < g{Uo), 

(1.6) 

^ < Po{Uo) < d,n{u,v) < h{Uo), 

where /3o, g and h are some constants depending only on Uq. Then there exists a characteristic speed 
a > large enough (see the condition \2. 5j) helow) such that the system il.l\) - [T7^) admits a unique 
solution {u^,v^) in C ([0, oo[, L;^q^(M)^) and there exists a constant C > which only depends on a and 
{uq,Vq), such that for any e > 0; 

\\v%t) ± Van=(t)||^ < C Vt > 0, 



(1.7) 



TV (u^it)) + TV (v^t)) < C Vt>0, 



\u%t + T) -u%t)\\i < Ct, Vt G 
\v%t + T)-v%t)h 



C 

< —T, Vt G 

e 



T G 



\\v^ {t + t) - {t)\\l < CuT, Vt>Z^,TGM+, 

where v > 0, and Cy only depends on a, {uq,Vq) and v. Moreover, there exists /3o > such that, 



(1.8) \\v^{tr)-A{u^{t,.))\\^ < e-- \\vl-A{ul)\\, + Ce. 

Finally, if Uq converges, as e goes to zero, to uq defined hy 111.5]) . then the sequence (u^,v^) converges to 
{u,A(u)) when e goes to 0, such that, for any v > 0: 

( ^ u C([0,oo);Li„,( 

(1.9) 

I V' A{u) C{[v,^)-Ll 
where u is the unique entropic solution to the Cauchy prohlem Ji.^l ). 

In this paper we propose a rigorous analysis of the Asymptotic Preserving scheme proposed by F. 
Filbet & S. Jin [23j and G. Dimarco &: L. Pareschi p/7j for a nonlinear relaxation. In other words, 
denoting by "P^ and respectively the relaxation and the equilibrium Cauchy problems, and "P^ and 
the corresponding discrete problems, where h represents the discretization parameter, independent of e, 
we will perform a precise analysis of the following Asymptotic Preserving diagram. 




The paper is organized as follows. We present in Section [2] an asymptotic preserving scheme for the 
relaxation model, and state the both convergence results of the Asymptotic Preserving scheme when the 
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relaxation parameter e goes to zero (Theorem I2.3P and next when the discretization parameter h goes 
to zero (Theorem 12. 4p . Then, we estabhsh some a priori estimates in L°° and BV on the numerical 
solution to the Asymptotic Preserving scheme in Section [3] in order to prove both the zero relaxation 
limit (Sectional) and the convergence of the scheme (Section [5]). Notice that these a priori estimates 
are uniform both with respect to the relaxation parameter e and the time step At and space step Ax. 
Finally, we present some numerical results in Section [6j 

2. Numerical schemes and main results 

We remind that when lZ{u,v) = v — A{u), where A £ C^{W,W) is a given function, the necessary and 
sufficient stability condition is given by the so called sub-characteristic condition [6\ f27j : 

(2.1) \A'{u)\ < V^. 

It means that the propagation speed of the equilibrium problem has to be bounded by the speeds of the 
relaxation system, which has to be dissipative. For more details about this case, we refer to |14 ^ 134 ^ 138]. 
Hence the sub-characteristic condition reads, in our case: 

duTZ {u,v) 
dvTl {u, v) 

For the sequel, we define for any > and q > 0, 




(2.2) 



' UiN,a):= + A, 
F{N,a):= sup \A{0\ 

m<U(N,a) 



^ V{N,a) :=U{N,a) + j=F{N,a) . 
We also denote by /(A, a) the compact set 

(2.3) I{N,a) := [-^/^V{N,a), ,/^V{N,a)]^ . 

Moreover, we assume that the initial conditions Uq,Vq are bounded independently of e in L°°(]R), such 
that: 

(2.4) Ao := max <! sup||uo||oo,sup||fol|oo > < oo . 

.e>0 e>0 



Consider any qq > and assume that the function TZeC^iMx M,M) satisfies ([L3|) and i^M)- We choose 
the characteristic speed ^/a > and the parameter /3 > such that 

g{V{No,ao)) 

V u ^ max \/ UQ , 

(2.5) 



I PoiV[No,ao)) } 
^ P = h{ViNo,ao)), 



where V is given by (j2.2p . 

Remark 2.1. Note that if we differentiate with respect to u the equilibrium equation 

n{u,A{u)) =0, 

we obtain 



(2.6) 



A'(u) 



dun{u,A{u)) 



d^n{u,A{u)) 

and thus recover the well known sub- characteristic condition in the case of semi-linear relaxation, namely: 

\A'{u)\ < 
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We present here the sphtting Asymptotic Preserving scheme and its relaxed version. We introduce a 
space time discretization based on a uniform grid of points {xj ^i/2)j £Z C M, with space step Ax, and 
discrete time t"" = nAt, n G N, for which the time step At satisfies the CFL condition: 

\/aAt 

(2.7) » < ^ - ^ 



< 1 



We denote hy h = {At, Ax) the discretization parameter. 



2.1. An Asymptotic Preserving scheme for the relaxation system. In this section, we design a 
numerical scheme for system (|l.ip - (ll.2p . by introducing a splitting between the linear transport part, 
and the nonlinear relaxation part, for which we will take advantage of the knowledge of the equilibrium 
p.3p . In the asymptotic regime e — )• 0, the differential equation (jl.ip becomes stiff and explicit schemes 
are subject to a stability constraints At = 0{e). Of course, implicit schemes allow larger time steps, but 
new difficulty arises in computing the numerical solution of a fully nonlinear problem at each time step. 
Here we want to combine both advantages of implicit and explicit schemes, that is, large time step for 
stiff problems and low computational complexity of the numerical solution at each time step. This is 
done, as said in the introduction, in the spirit of Asymptotic Preserving schemes for nonlinear relaxation 
problems introduced by F. Filbet & S. Jin [23] and G. Dimarco &: L. Pareschi |17] . 
Thus we construct a numerical solution {uf^,vf^) to ()l.ip - ()1.2p in x M as follows 



(2. 



neN ji 



nSN ji 



where Cj =]xj_i/2jXj^i/2[ are the space cells and the sequences {u'j)nj and {Vj')nj depend on e and are 
given below. 

First, initial data are computed as the averaged values of (jl.2p through each space cell: for all j in Z, 



(2.9) 



1 

Ax 



Uq{x) dx , 



-L X/o(-)dx. 



Then, in order to discretize the system (jl.ip - (jl.2p . we apply a splitting strategy into a linear transport 
part and a stiff ordinary differential part as follows. The first part consists to apply a time explicit scheme 
combined with a finite volume method to the following linear differential system 

{dtu + dxV = 0, 
dtv + a dxU = , 

and then the second part deals with the stiff ordinary differential equations 

dtu = 0, 

(2.11) , 1 

dtv = — Tl{u, v) . 



We first approximate the linear transport part (12.10p . that is, for a given (u" , f"), we compute 
time t"'~^^ with a standard Finite Volume scheme, that is, for all j E Z, 



(2.12) 



3 

n+1/2 



At a DhUj , 
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where D^Vj' and aD^vJ- are discrete derivatives with respect to x of f and u, given for instance by a 
Lax-Friedrichs scheme, namely: 

1 



2 Ax 



2A^ L« l^i+i 



Remark 2.2. Of course, there is a wide range of possible choices for the numerical fluxes. As we will 
see below, the main property of the numerical scheme for the linear transport term that we require is the 
TVD (^Total Variation Diminishing^ property, namely, for all n E N, 



[ TViv 



n+l/2\ 



< TV{v'- 



where TV{u) :- 



Hence, the second part of the sphtting only consists in approximating the nonlinear ordinary differential 



equation (I2.1ip . for all j G Z, starting from {u 



n+l/2 n+l/2-. 



To this aim, we use the decomposition 



n{u, v) = [n{u, v)- f3 {v- A{u))] + P {v- A{u)) , 
where /3 > is a parameter such that 

< sup dvTZ{u, v) < (3. 

{u,v) 

Then, we apply a time exponential scheme on the dissipative part and get the following numerical scheme: 



(2.13) 



,n+l 



U 



V 



n+l/2 

'3 

n+l/2 
j 



n+l/2 



u 



n+l/2 n+l/2 
1} ■ 



This numerical scheme (|2.12p - (j2.13p allows to define a sequence 



eZxP 



2.2. Convergence results. We first establish a convergence result on the asymptotic behavior of the 
numerical solution to (j2.12p - (j2.13p when e tends to zero and h = (At, Ax) is fixed. 

Theorem 2.3. Assume that the initial conditions {uq,Vq) are bounded independently of £ in BV{E.) and 
such that the assumption {2.4\) is satisfied. Consider IZ E C^(M x ]R,M), which satisfies hl.3\) - [LBjl and 
the characteristic speed \fa > and the parameter /3 > are given by \2. 5|) . Then, the solution {uf^,vf^) 
given by i2. 8\) to the scheme i2.12\) - r2.13\) with the initial data 112. 9\) . converges in -L^(M), as e — )• 0, to a 
numerical solution {uh,Vh), that is, 

\\uUt)-Uh{t)\U + \\vm-Vhmi < Cte-^"^'/' + 
where {uh,Vh) is a consistent approximation to the conservation laws \1.4^ with = A{uh), 

Uh{t,x) := ^ 1C,(X) l[in_in + l[(t), 

and the sequence 

(2.14) u]+^ = + AtDhA{u]), i e Z, n > 1, 
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with the initial data 

(2.15) u^j = ^ Mx) dx, 

where uq is given by U.5\) and 6^ = Uq — A{vq). 

The proof of this theorem corresponds is analogous to the one corresponding to the continuous problem. 
The fact that we are able to establish uniform BV bounds on the numerical solution allow to get error 
estimates with respect to e. 

On the other hand, we propose a convergence result of the Asymptotic Preserving scheme when 
h = (At, Ax) goes to zero and e is fixed: 

Theorem 2.4. Assume that the initial conditions Uq,Vq are bounded independently of £ in BV{M) and 
such that the assumption ^2.4\ ) is satisfied. Consider TZ G C^(MxM,M), which satisfies il.3\ )- [T7S\) and the 
characteristic speed y/a > and the parameter /? > are given by \2. 5|) . Then, the solution (u|,f|) given 
by / fO|) to the scheme \2.12\) - [KW\) and the initial data IO|) . converges in L]^^ (M+ x M), as h ^ (0, 0), 
to a weak solution {u^,v^) to the relaxation Cauchy problem hl.l\) - [L^) . More precisely, we have 




\\ul(t) - n^(t)||i + \\vm - v^{t)\\, < - ( At ( + 1 1 + Axi/2 

where 5^ = Uq — A{vq). 



3. A priori estimates 

We first define precisely the parameter /3 and the assumption on the characteristic speed to ensure 
the stability of the scheme ()2.12p - ()2.13p and prove estimates on the solution to the relaxation problem 
which are uniform with respect to e. In following section, we drop the subscripts e for sake of clarity and 
investigate the stability property of the Asymptotic Preserving scheme ()2.12p - ()2.13p . 

3.1. A priori estimate on the relaxation operator. Let us focus on the second part of the scheme 
devoted to the approximation of the relaxation source term and give a technical lemma, which establishes 
a quasi-monotonicity property on the operator G^^s with s > 0. In order to do this, we will rather consider 
the equivalent formulation on the diagonal variables w and z. Let us rewrite the splitting scheme on 
these variables. For u and v given, we set 



w = —V — Wau. 



(3.1) 



+v — -y/a 



au . 



Therefore, the linear transport scheme (|2.12p written for (w, z) exactly coincides with an upwind finite 
volume method: for all j G Z, 



(3.2) 



J ^ ^ Ax 

whereas the nonlinear stiff part (|2.13p yields, for all j G Z 



(3.3) 



n+l n+1/2 , ^ ( n+l/2 n+1/2 
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with 



z — w , I w + z 
A 



e 



2^ 

w + z z — w 
' 2^ ' 2 



The main result of this section is to estabhsh the quasi- monotonicity of the operator Gg^s^ which wiU 
lead to L°° an BV estimates. 

Lemma 3.1. Assume the function TZ G C^(M x R,M) satisfies il.3]) - [T^) and choose a > 0, (3 > such 
that h2. 5|) is verified. Then, 

(i) the suh- characteristic condition is satisfied for all {w,z) £ I{NQ,aQ), that is, 

dun, 



(3.4) 



a,7^ 



{u,v) 



where := —{w + z) and 2v := z — w ; 

(ii) for all £, s > 0, the source term operator G^^s is quasi-monotone on the compact set I{No,ao), 
that is, 



(3.5) 



T < dujGe,siw,z) < 0, y{w,z) G /(A^o,ao), 



< d,Ge,s{w,z) < 1, y{w,z) e /(iVo,ao); 

{Hi) consider for i = 1,2, {w'j^'^^ , z^^^) two solutions to \3. 3\) corresponding to two initial data 

{u\~^^^'^ ,z^^^^'^) G I{NQ,ao). Then there exist w and z G M such that \w\, \z\ < ^/aV{NQ,ao) 
and 



(3.6) 



W 



n+1/2 n+l/2\ 



Wr. 



2 ' j {1 + dyjGe,s{w,Z^ 



n+1/2. 



+ U 



n+1/2 n+1/2 



'1 



z'-T^'"'] dzGe^siwl'^^^'^ ,Z) 



n+1 _ n+1 



n+1/2 n+1/2 
'1 ^2 



— W 



n+1/2 n+1/2 



Wo 



1 - dzGs,s{w2~^^^'^ , z] 



dyjGe,s{w,z'^^^^'^) . 



Proof. For any A'^o > and > 0) we first observe that for (w, z) G I{Nq, oq), 

II |t« + 2:| 



2VE 



<V{No,ao). 



Therefore, using the assumption (II. 6p and the definition (12. 2p . we get that 



{u,v) 



g{V{No,ao)) 
< ^ r- < va. 



^o{V{No,ao)) 

which proves the first assertion (i). 

Now we prove (ii) the quasi- monotonicity property of Gs^s- Computing the partial derivatives of Gg 
it yields for all s > 0, 



2 \ 



dwGs,. 



1 _ ( 1 + ^ ) e-"'" 



A.-...(^^,,,,), 



^ 2e 



-dull 
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Hence, from Remark l2.1l and the sub-characteristic condition ()3.4p . we obtain that for all {u, v) G I{Nq, oq) 

dujGe,si'w,z) < and dzGa,s{w,z) > 0. 
Moreover, still using condition (j3.4p . we also get for all (wjz) £ I{NQ,aQ) 



dwGe,siw,z) > 



dzGe,s{w,z) < 



1 



/3s 



+ dvTZiu, v) - e 



-I3s/e 



Now since \u\ < V{No, oq) and from the choice of the parameter /3 in (|2.5p . it yields that \dvTl{u, v)\ < 
p. Therefore, we conclude that 

-1 < dwGe,siw,z), dzGe,s{w,z) < 1, ^ {w , z) £ I {Nq, ao) ■ 

Finally (iii) follows from a first order Taylor expansion of G^^s- D 

This Lemma allows to obtain the following comparison principle. 

Corollary 3.2. Consider for i = 1,2, two initial data {w'^~^^^'^ , z^^^^"^) G I{NQ,ao) satisfying the mono- 
tonicity condition 



w 



n+l/2 
1 



< W. 



n+l/2 



n+l/2 n+l/2 



and Zi < z. 



Then, the numerical solution {w^ , ), given by i3. 3\) corresponding to the initial data (w. 
for i = 1,2, satisfies 



n+l/2 n+l/2 



W 



n+1 



< W. 



n+1 



and z 



n+1 
1 



< Z. 



n+1 



Proof. Starting from the equality (|3.6p . it yields to the result applying the estimates (jS.Sp . □ 

3.2. estimates. We establish a uniform bound on the numerical solution to the scheme (j2.12p - (j2.13l 
with respect to the time-space step h = (At, Ax) such that (j2.7p is satisfied. 

Proposition 3.3. Consider any ag > and 

Nq = max < sup||iio||oo, sup||wo||( 

U>0 e>0 

We assume that the function TZ £ C^(M x M,M) satisfies lll.3\) - [T^} and choose a > 0, /3 > such that 
f2.5\} is verified. Then, for all n £ 'N 

||n"||oo < V{No,ao), ||i;"||oo < V^V{No,ao). 

Proof. We proceed in two steps and first construct a particular solution {w^,'z^) G to the scheme 
(|3.2p - ()3.3p which is uniformly bounded, then we apply the comparison principle on the compact set 
I{NQ,ao) to prove an bound on (u",f"). 

To this aim we choose Rq = {1 + ^/a) Nq and construct a numerical solution {w^,'z^) to ()3.2p - ()3.3p 
corresponding to the initial data (vf,!^) = {Ro,Ro), which does not depend on the space variable so 
that the transport step (j3.2p is invariant. Then we apply the relaxation scheme (|3.3p . which yields 



w 



where (u^,Jf^) are only given by 



13 

l + - I e 



0At/e-n~l ^ fi _ fi + ^] e-^^'A A (nO 
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Then, we proceed by induction to show that: 

Vn G {0,...,iV}, (lI7",z")G/(iVo,ao). 

We assume that (uJ""-*^, z""-*^) G I{NQ,aQ), for some n > 1. Let us prove that (itJ",z") G I{NQ,ao). On 
the one hand since = TiP, it yields 

1 



^ oo — oo _ 



1 + 



A^o < U{No,ao), 



where U{No,ao) is given by (12. 2p . 

On the other hand, using a first order Taylor expansion of the source term TZiiiP, .), we get that there 
exists v^~^ G M such that 



e 

Therefore, denoting by G M, the real number such that 



A, := ( 1 + ^-Wn',v') ^ ^-M./., vfc G N, 



with \v^\ < F(iVo,ao), hence we get 

v'^ = Xn-iv""-' + (l-A„_i)A(nO) 

and since iP = 



n-l \ 

1 - n A J A{u") 

k=0 J 



V = I 

Moreover, using that \v°\ < U{No,ao) and < ^/aV{No,ao), for ah G N, we get from (ll.6p and ()2.5p . 

< (^1 + At) e-/^^*/^ < Afc < (^1 + ^ ) e-'^^*/^ < 1, Vfc G N. 

Therefore, ||l'"||oo < F{No,ao) and 

lltx^^lloo, ||^"||oo < F{No,ao) + {l + ./^)No<V^V{No,ao), 

that is, G /(iVo,ao)- 

Furthermore, starting from the following initial datum {wP,z^) = {—Ro,—Ro), we can also construct 
another particular solution {uf^^z^) G I{Nq, gq) for all n G {0, . . . , A^}. 

Now, we proceed to the second step which consists in applying the comparison principle of Corollary 
13.21 to prove an estimate for any initial data u^, G -L°°(M) given by ()2.9p . From the definition of 
A^O) we have 

||u°||oo, |b°||oo < iVo. 

Then, we have for the initial data {w^, z^) 

lk°||oo, Ikloo < {l + V^)No = Ro <V^V{No,ao). 
In other words, we have initially: 

Thus, we proceed by induction and assume that 

u;" < li;" < IZJ", < z'^ < z". 

We first consider the linear transport step (13. 2p to {w^'',z"') and get that 
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Thus, we apply Corollary l3.2l to the two solutions to (I3.3p associated to the initial conditions (w"^^^^, z"^^^^) 
z""'"-'^/^) and {w2^^^'^ , Z2~^^^'^) = (and then (l(J",z")), we have 

which finally gives for all n E N, that {w"',z"') G I{NQ,aQ). By construction of {vP,v^) we have proved 
that 

h"||oo < V^(iVo,ao), lb"||oo < ^V{No,ao). 

□ 

3.3. BV estimates. In this section, we obtain a BV estimate on the numerical solution to the scheme 
(j2.12p - (j2.13p with the time-space step h = (At, Ax) such that (j2.7p is satisfied. 

Proposition 3.4. Assume that uq^vq are uniformly hounded with respect to e in BV{M.). For any oq > 
and Nq = max < sup||no||oo) sup||wo||oo \, we assume that the function TZ G C^(M x R,M) satisfies il.3\) - 

U>0 £>0 J 

lll.6\) and choose a > 0, /3 > such that /i2. 5|) is verified. Then, for all n £ 'N, we have: 

ry(w;"+i) + rF(z"+i) < TViw"") + ry(z"). 

Proof. First we note that u^, £ BV{M.), then by construction, w^, z^' G BV{M.). To prove the BV 
estimate, we proceed in two steps. On the one hand, using the TVD property of the upwind scheme, we 
get that 

rF(ii;"+i/2) < Ty(u;") and TV{z''+^l'^) < Ty(z"). 
On the other hand, we apply the nonlinear relaxation step (13. 3[) and from Lemma l3.ll iiii) with w"^^^^ = 

n+l/2 n+l/2 n+1/2 , n+1/2 n+1/2 n+1/2 n+1/2 • , , r ■ ^ n? 

Wj , Zi = z- and W2 = ■, '^j+i ' yields for any j G i^, 

I n+l n+l| I I n+1 n+li ^ 1 n+1/2 71+1/2, , , 71+1/2 n+1/2 , 

Summing over j G Z, we get that 

ry(u;"+i) + ry(z"+i) < tv{w''+^/'^) + tv{z''+^/'^) 

□ 

4. Trend to equilibrium (proof of Theorem I2.3p 

For a sequence u = {uj)j,^z we set 

||u||i := Ax \uj\. 

In this section we first focus on the asymptotic behavior of the numerical solution to (j2.12p - (j2.13p when 
e goes to zero or when times goes to infinity. Then, we prove that this numerical solution converges 
to a consistent approximation of the conservation laws (|1.4p when e goes to zero. Roughly speaking, it 
corresponds to the limit V\ — )• V\., when e — )• 0. 

4.1. Asymptotic behavior. In this subsection, we drop the subscripts e for sake of clarity and estimate 
the deviation to the local equilibrium = v'^ — ^(u"). 

Proposition 4.1. Assume that uo,vq are uniformly bounded with respect to e in BV{W.). For any oq > 
and Nq as before, we assume that the function TZ G C^{W x M, M) satisfies lll.3\} - p^) and choose a > 0, 
(3 > such that Ii2. 5\) is verified. Then the deviation from the equilibrium, 6 = v — A{u) satisfies for all 
n G N and all e > 

||^n+l/2||^ < + CAt, 



(4.1) 



< e-^o*"/^ \\6^\\^ + Ce. 
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where C > is a constant only depending on the parameters a, (3q and the BV norm of the initial data. 
Moreover, if £ < At then we get 

(4.2) < e-^"*"/^ \\6^\\^ + CaAte-l^°^*/^. 

Proof. We set for j G Z, n G N the sequence of deviations from the equilibrium: 

We first consider the transport step (j2.12p of the numerical scheme: for all j £ Z, we apply a Taylor 
expansion to A, then there exists ^" such that |^"| < V{NQ,ao) and 



At 

2 Ax 



Thanks to the uniform BV estimate, proved in Proposition 13.4^ the sub-characteristic condition 

and the TVD property of the numerical fluxes we get the first estimate (|4.ip . by multiplying by Ax and 
summing over j €z Z: 

(4.3) < n + AtCa [TV{v°) + ,/^TV{u°)] , 

where Ca > is a constant only depending on a. 

Then, we consider the second step of the numerical scheme (j2.13p . On the one hand, since n"^^ = 
it yields 



en+l _ rn+1/2 
j j 



1 + 13 



At 



At 



On the other hand, applying a Taylor expansion, since Tl{u, A(u)) = we get that there exists r/ such 
that \r]\ < -/aV{No,ao) and: 



7^(n;+^/^^;;+^/') = a.7^(n;+^/^ 7?) 5;+^/^ 



Hence, we have 



^n+l _ j"+l/2 



Therefore under the assumption (11. 6p . we set for all s > 

for which we easily show that for all s G we have that e^^ < g {s) < e~^°^^^ . Hence, for s = At/e 
(4.4) P'^'Wi < e-^«^*/n|5"+i/2||^_ 

Finally, gathering (j4.3p and (j4.4p . we obtain that there exists a constant Ci > depending only on a, 
ry('u°) and TV{v^) such that 

By induction, we easily get 

g-/3o At/e 



(4-5) ll^lli < e-^«*"/^ 11^1, + C.At-_^_^^^^/^. 

To conclude we only observe that xe~^ < 1 — e~^, for any x > 0, then it gives the second estimate of 
(|4.ip : there exists a constant C > 0, only depending on a, Po, TV{u^) and Ty(t''^) such that 

linii < e-**"/^ \\6'\\^ + Ce. 



AN ASYMPTOTIC PRESERVING SCHEME FOR RELAXATION SYSTEMS 



13 



Moreover, when e < At, we again start from the estimate (j4.5p and note that 1/(1 — e~'^°^*/'^) < 
1/(1 — e~^°). Thus, there exists another constant C > 0, only depending on a, /3o, TV{vP) and TV{v^) 
such that 

which gives (j4.2p . □ 



4.2. Proof of Theorem 12.31 We are now ready to perform the asymptotic analysis of the numerical 
scheme (|2.12p - (j2.13p when e goes to zero. 

Let us consider the numerical solution {uf^,vf^) to the scheme (j2.12p - (|2.13p written in the form (j3.2p - 
(I33D with 



'au 



4 = +'"h 



h ' 



aul, 



such that 



' ■wUt,x) = ^^u;^'"lc,(x)l[tn^in+i[(t), 



new J 6 



4( 



Let us also define {w^, z^) the numerical solution to the scheme (|3.2p - (|3.3p in the asymptotic limit e = 0. 
Wh and Zh are given by 

Wh{t,x) = ^ ^u;" lc, (x) l[jn^t„+i[(t), 

= ^^zjlc,(x)l[tn^tn+i[(t), 

neN jez 

where (w", z") is given by (j2.14p . To obtain an error estimates we rewrite the values [w^ , z'j)(^n,j)eNxZ as 
a perturbation of the numerical solution to ()3.2p - ()3.3p with a fixed value of e > 0, 



w 



ra+l/2 



n+1/2 



Ax 
r- 



and then 
(4.6) 



w 



n+l 



n+1/2 



n+1/2 n+l/2\ 



) - At £:;(£), 



^n+l _ n+1/2 



n+1/2 ^n+1/2 



+ Atf;^(e). 

where AtiS"(e) represents the consistency error of the operator Ge^At with respect to e, that is, 



n+1/2 n+l/2\ 



n+1/2 n+l/2\ 



Therefore, we apply Lemma [3TT] (ii) and (iii) in (j4.6p . with (t(;i,zi) = (t(;|,z|) and {w2,Z2) = {wj,Zj), it 
yields 

I e,n+l n+li , i e,n+l n+li ^ i e,n+l/2 n+l/2| , , e,n+l/2 n+l/2| 

\W^' -W^+\ + \z^' -Z + \ < ' -W^ + '-^3 M 

+ 2|At£:;(e)|, 

and by linearity of the transport scheme ()3.2p . we have for all n > 



I e,n+l n+li , i e,n+l n+li ^ i e.n ni , i e.n ni 



+ 2 |At£:;(e)|. 
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Thus, multiplying by Ax, summing over j £ 1^ and applying a straightforward induction, we get the 
following stability result 



^ Ax (\w'.''' - <| + - z^\) < ^ Ax (\wf - + {z 



,s,0 _ 



n-l 



+ At Ax £f{e) 

k=ojez 

It now remains to evaluate the error <f"(e). Using that for any {u,v) G I{NQ,aQ), the function TZ G 
CHM+,M) such that /3o < d^n{u,v) < (5, we have 



At \8^{e) 



-/3Ai/e 



< g-/3oAt/£ 



n+1/2 
3 



n+1/2 n+1/2 - 

- \(u^ 

n+l/2> 



V. ■ -A(y:'^''^)\ 1 + 



/3At 



At 



+ —71 



n+1/2 n+1/2 



3 



3 



Thanks to the estimates (j4.ip and (]4.2p in Proposition 14. 1 1 on the deviation applied to ti^+^Z^ — 
which is also valid in the asymptotic e — )■ 0, it yields 

' + C At if n = 0. 



|^-+l/2_^(^n+l/2)||^ < ^ 



CAt 



if n > 0. 



Then, we get for A: > and e < At, 

Y^AxAt £^{e) < e"*^*/" + C At] 



Hence summing over < A: < n, it gives 



< e-^°^*/^ [\\6\ + Ct^]. 



^J^AxAt £^^^, 
k=o jez 

Finally, we get the estimate 

- wunwi + \\zun - zkinwi < - w.mh + \\zm - zmwi 

and the result follows {uf^,vf^) — )■ {uh,Vh), when e goes to zero. The proof of Theorem 12 .31 is now complete. 

5. Proof of Theorem 12.41 

In this section, we prove the convergence of the relaxation Asymptotic Preserving scheme. As in the 
stability analysis of the relaxation scheme, we will rather consider the diagonal variables w and z and 
drop the subscripts £ for sake of clarity when it is not necessary. 

5.1. Consistency error. Consider {w,z) the exact solution to (|l.ip - (|1.2p with (|3.ip . Unfortunately, 
this solution is not smooth enough to study the consistency error, then we introduce a regularization 
{ws,zs) given by 

ws{t,x) = w-k ps{t,x), 

zs{t,x) = z-k ps{t,x), 
where ★ denotes the convolution product with respect to x G M and 

ps{x) = ^p(^) and p G C^{R), p > 0, j p{z)dz = 1. 
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(5.1) 



dtws + \/adxWs = H — TZs{u,v), 



dtzs - VadxZs = — TZs{u,v), 



with TZs = TZ-k ps and {u,v) solution to ()l.ip - ()1.2p . Therefore, applying the Duhamel's representation, 
the solution can be written as 



(5.2) 



( 1 Z"^* 

u;5(r+\x) = ws{t",x - ^/aAt) + - ns{u,v){e + s,x - y/a{At - s)) dt, 

^ Jo 
1 f^* 

zs{t"-^^,x) = zs{t'^,x + ^/aAt) / ns{u,v){t''^ + s,x + ^/a{At - s)) dt. 

V e Jo 



Then we set 
(5.3) 



^ Ax 



ws{t'^,x)dx, z"" = — 
■' Ax 



zs{t"-,x) dx. 



Integrating (j5.2p over x £ Cj and dividing by Ax, it yields 
(5.4) 



J J 

n+l _ ~n+l/2 
^3 ~ ^3 



~n+l/2 ~n+l/2 



G/ -11-1-1/ ^ ~r 



Ateij + At 82 J, 



with 



-n+l/2 ~n r~ I ~n ~n \ 



n+1/2 _ „ 



Ax 
At 



The consistency errors related to the transport operator £^ are respectively defined by 



A.pn _ ^l,j+l/2 ^l,i-l/2 . „ 



^3j+l/2 ^3,i-l/2 



Ax ' •^'^ Ax 

where j_|_i/2 ^3 j+i/2 ^'■^ consistency errors of the numerical flux and are given by 



(5.5) 



£ 



i,i+i/2 



/ ^^(t", Xj_|_i/2 — •s)ds + \/aAtw^. 
Jo 



^3j+i/2 = + / ^5(*">3;j+i/2 + s)ds - ^/a At 
whereas the consistency errors At j and At £2j correspond to the stiff source term and are given by 

At£:£,. = +^ / -TZs (n, v) (t" + s,x- ^^(At - s)) ds - Ge,At (w^'-'^^ z']^^''') dx. 



n+1/2 ^-n+l/2\ 



Atf4"' . = --^ / / -7^^ (u, w) (t" + s, X + V^(At - s))ds - At (w 
Ax Jc^ Jo £ V ■ 

We then evaluate successively each consistency error term. On the one hand, we prove the following 
consistency error for smooth solutions, which is related to the transport approximation. 
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Proposition 5.1. Let {w,z) be given by i3.1\} . where {u,v) is the exact solution to M.l\j - [L^j and such 
that w, z € L°° , BV (M)) . Then the consistency error related to the transport part satisfies 



Y,Ax [\£^J + \£l^\] < C— {TV{w{t-)) + TV{z{t^)). 

Proof. We first study tlie consistency error for w G i?y(M)). We perform a simple cliange of 

variable in (15. 5p . which yields since ^/a/S.t = A Ax, 



e 



-X ws{t'^,Xj^i/2- Xs)ds + X ws{t"-,Xj+i/2 - s)ds, 
Jo Jo 



Ax rs 



10 JXs 

Therefore, since ws is smooth we have 

f Ax rs 

. A / / / I 



A/ / dxWs{t'^,Xj+i/2 - r)drds. 



JXs 



i-3/2 



dxWs{t''^,Xj+i/2 -r) - dxWs{t",Xj_i/2 - r)drds 
d'^^ws{t^,x)\ dx. 



By multiplying by Ax and summing over j G Z, we get an estimate for a smooth solution ws{t^) G M^^'^ 

J] Ax \£l^\ < 2V^Ax\\dl,ws{tn\\ 
jez 



1- 



To achieve the proof, we need to estimate with respect to w and ps- Using the convolution 

properties, we easily get 



which allows to conclude that 



J^Ax < C^TV{w{f^)). 



Using a similar technique, we also get for a smooth solution z G L°°(M"'", BV{M)), 

Y.Ax\£l^\ < C^TV{z{t-)). 



□ 



On the other hand, we treat the consistency errors £2^ and £2j, which are related to the stiff source 
term. 

Proposition 5.2. Let {uj,z) be given by (3. where {u,v) is the exact solution to ( li. jji -f TOj) . Assume 
that w, z £ L°°(R+, i?y(M)). Then there exists a constant C > 0, only depending on u and v such that 
the consistency error related to the stiff source part satisfies 



and 



J2Ax\£SJ < C 



^Ax \£lj\ < C 



At (^.^.t^/An^A ^ Ax 



6 

+ - 

£ e 



At 



15° 



Ax 6 



e £ 
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17 



ul; = ^ — -pr^ and ti" = — -. 

^ 2Va ^ 2 



Therefore, we spht the consistency error ifj j as 

^2 J — ^21 J "T ^22 J "T ^23 J "T ^24,j ^25,j' 

with 



At^:2"2,i 



e 



~n+l/2 . /~n+l/2 



1 r z"^* 

At£^^j = —— / 7e5(u,t;)(r + s,x- Va(A^-s))-7^^(n,^>)(^,x- Va(At))dsd2;, 
/ 7^^(n, v){e,x - Va{At)) - TZ{u, v){t'^, x - Va(At))dx, 



At 

eAx 

At 



^^'^ eAx Jc^ 

On the one hand, the two terms £21 j and 1^22 j can be easily evaluated using a Taylor expansion of 
s ^ e~^^^^, it yields 



At\£^,J < 



~n+l/2 . [ ~n+l/2 



Using that 7^(^i,^(n)) = and 7^ G C^(M^]R) with d^Tliu.v) < (3, we also obtain that 



Therefore, from (j4.ip in Proposition 14. H we have 



~n+l/2 , /~n+l/2 



(5.6) 



5:Ax[|£:2\,| + |f2\,l] < C^(e- 



On the other hand, we proceed to the evaluation of the terms £23 ^mj ^^^^ ^25 j- First, for s € [0, At], 
we set 

^5,x{s) = [Tl{u,v) -k ps] (t" + s,x- ^/a{At - s)). 

Then, from (jl.6p and (j2.5p . we know that \du'R.{u, v)\ < y/a (3 and \dyTZ{u, v)\ < /3, for any {u, v) G oq), 
we obtain 



■^n, ' ^ JR Jo Jo 



dx, 



At r r*"^' 

< C— / / {\dtus\ + \d^us\){t,x)dtdx 



+ c 



At 



i\dtvs\ + |«9x?^5|) it,x)dtdx. 
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Thus we can use the estimates on the continuous relaxation system hsted in Theorem ll.il Indeed, since 

dtus = -d^vs, 

dtvs = -adxus - ^TZs{u,v), 
we obtain, by applying a first order Taylor expansion of TZ, the inequalities 

i\dtus\ + \d,us\){t,x)dx < TV{u{t)) + TV{v{t)), 

{\dtV5\ + \dM){t.^)<^x < C (TV{u{t)) + - \\{v - A{u)){t)\\^ . 
Hence, integrating over t G and using (|1.7p and (jl.Sp . we get: 



(5.7) J^Ax <C—\ TViuin) + TVivin) + 

where C > only depends on ^/a and /3. 



-/3ot"/£ 



Now we treat the term £"24 j using the smoothness properties (|1.6p and (j2.5p on TZ, it gives 
VAx|f2\jl = - / [[J^{u,v){e,x-y-^/^At)-niu,v){e,x-V^At)]p5{y)dy 



dx, 



< ^ 



|n(r,x)-n(r,x-2/)| + \vie ,x) - v{e , x - y)\] ps{y) dy dx. 



Thus, applying Fubini's theorem the BV estimate on the exact solution ()1.7p and the value of the integral 
of PS , we get 

6 

e 



(5.8) 



5^Ax|£:2\,| < C-[TV{u{e))+TV{v{n)]. 

Finally, to deal with the last term ■ , we split it in two parts 

VAx|£'2"5| < - 1 \n{u,v){e,x-y/^M)-n{u5,V5){f',x-y/aAt)\dx 

~7i ' £ JR 



+ -J2 I '^ius,Vs){t'^,X - ^/^At) -n(u 



n+1/2 ~n4-l/2 

3 



dx 



and treat the different terms as for <?24 j' S^^ ^^'^ ^'^^^ 

|7^(n,^;)(^",x) -7^(n^,^.^)(^",x)| dx < C6 [TV{u{e)) + Ty(y(t"))]. 
and for the latter one using the BV estimate on the exact solution (jl.7p . 



V / 7^ (U5, i;^) (t", X - VaAt) - 7^ 



~n+l/2 -n+1/2 



dx 



Thus, we have 
(5.9) 



< cAx [Wd^usinwi + \\dxvs{n\\i,]. 

^Ax |f2\,| <c(^- + ^\ [TV{u{n) + TV(v{n)[ 
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Gathering (j5.6p . (j5.7p , (j5.8p and (j5.9p . and finally using the uniform bound on the BV norms of {u,v) 
given in ()1.7p . it yields 



At 



-Pot 



6 

+ 11+ — + - 

e lee 



Using the same arguments we also prove that 



e I e e 



□ 



5.2. Convergence proof. Now we perform a rigorous analysis of the numerical scheme (|2.12p - ()2.13p 
when h = (At, Ax) goes to zero. We consider the numerical solution {uj^,vD to the scheme ()2.12p - (j2.13p 
and the exact solution to (jl.ip - (jl.2p and define {w^^z^) using ()3.ip . Then we denote by 

u,n = ^ f w%e,x)dx, z'^ = ^ f z%e,x)dx 
^ Ax Jc^ ^ Ax Jc^ 



and (tLi", 2")(j_„)g^xN ttis numerical solution given by 



Thus, 



Ax [ \w'; - w]\ + - zjl ] < ^ Ax [ \w] -w';\ + \z^ - ~z]\ 



where {w'j,z^)Q^n)eZx'N is given by (|5.3p . On the one hand, we estimate the second terms of the right 
hand side using the convolution properties and have 



(5.10) 



YAx[ 



- w]\ + \z^ - z]\ ] < C5 [TV{u) + TV{v)]. 



On the other hand, we apply the consistency error analysis to estimate the first term of the right hand 
side. Applying (j3.5p - (|3.6p established in Lemma [3TT] with (wj^Zj) and (wj^zj), it yields 



'^AxlwJ 



^ A I ~n+l/2 n+1/2. ( . , p. r< ( n+l/2^ 

< 2^Ax\w^ ^^j I [^ + 9wGe,At{Wj,Zj 



+ '^AxlzJ'^^^'^ - Zj'^^^'^ldzGg^Atiw, ' ' ,z 



n+l/2 



and 



YAx\zJ+'^ - zJ+'^\ < J^Ax|5j 



~n+l/2 n+l/2| 



3 



1 - d;,Gs^Atiw]^^^^,Zj 



n+l/2 n+1/2 



n+1/2-. 



+ ^AxAt[|f3^^.| + |fr, 
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Summing the two inequalities and using that the scheme ()3.3p is TVD, we get the following inequality 



je2 



Therefore, 



^Ax 



<Y^Ax[\S^-^\ + \u/^-v.^\ 



3& 



A:=0 j& 

Finally the consistency error analysis performed in Propositions 15.11 and 15.21 we have taking 5 = Ax^/^ 



r^i + i*r^-<^i] <ea- + 



(5.11) 



C 



+- At(At + e) 



15° I 



+ 1] +f- 



Ax + eAx^/2 + ^^1/2 



Gathering (|5.10p and (|5.1ip . the right hand side converges to zero when h = (At, Ax) goes to zero, which 
proves the convergence of the numerical solution (I2.12p - (l2.13p to the exact solution to (ll.ip - ()1.2p . 

Therefore, from the smoothness of the exact solution and the initial data {uq^vq), it proves that for 
any discretization parameter /i, for all e > 0: 

c ( /||<5or 

K(t, x) - u%t, x)\ + \vl{t, x) - v%t, x)\dx < — I At - — - 



I 



+ 1 + AxV2 



6. Numerical simulations 



6.1. Generalized Jin& Xin model. We consider as a first numerical test the Jin & Xin model with a 
nonlinear source term 

dtu + dxV = 0, 



(6.1) 

with Tl{u, v) given by 



dtv + adxU 



TZ{u, v) 



1 



n{u,v), 



V — u 
V? + v"^ 



We compute an approximation of the error for different meshes in space and time in order to evaluate 
the order of accuracy of the numerical scheme and different regimes corresponding to small and large 
values of the relaxation parameter e > 0. Therefore we choose two different initial data : the first one is 
a smooth initial datum given by 

(6.2) no(x) = sin(27rx), vo{x) = Q, Vx € (0, 1) 

whereas the second one is discontinuous 

' 0.5 if X < 0, 
0.125 ifx>0, 



(6.3) 



Uq{x) 



and vq = 0. 

We define an estimate of the numerical error by 



£i{h) = 114 -n|;,||i + \\vl-vlh\ 



1) 
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where h = {Ax, At) is the discretization parameter and || • ||i is the discrete norm. In Figure [H we 
present the curves corresponding to the order of accuracy with respect to h, computed from £i{h) and 
observe that the order is relatively closed to one when the solution is smooth, whereas it decreases to 
1/2 when the solution is discontinuous which occurs when e goes to zero or when the initial datum is 
discontinuous. Furthermore, as it has been proved in this paper, the numerical error is not too much 
affected by the variations the relaxation parameter e. 




Figure 1. Test 1 : error (a) smooth initial data i6.2\) and (h) discontinuous initial 
data for different regimes with initial data lid. 8\) . 



6.2. The Broadwell model. We now apply our scheme to the Broadwell model, which can be seen as 
a simple one-dimensional lattice Boltzmann equation, with only three discrete velocities [HI [33]. The gas 
is defined by a density function in phase space satisfying the equation 



(6.4) 



(/o - /+/-), 



dtfo — ^ (/+ /- - /o ) , 



(/o -/+/-)• 



Here /+, /_ and /o denote the particle density distribution at time t, position x with velocity 1, —1 and 
respectively; whereas e is the mean free path. 

We can rewrite the system with fluid dynamical variables. First, we define the density p, the momentum 
m and z as 



(6.5) 

and the system (16. 4p reads as 
(6.6) 



P -- 
m 

z - 



dtp + d^m 
dtm + dxZ 

dtz + dxTn 



/+ + 2/o + /_, 
-- f+- f~, 
/+ + /-• 



0, 
0, 



m 
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The local equilibrium corresponding to ()6.6p is then defined hy z = A {p,m), where 

,2^ 



A{p,m) 



P + 



m 



P 



Hence, when e goes to zero, we obtain the following system of equations 

dtp + d^m = 0, 

(6.7) 

dtm + dxA{p,m) = 0. 

We still perform several numerical simulations on this model in order to verify the order of accuracy and 
the asymptotic behavior of the numerical solution for different values of e > 0. 

Test 1 : Approximation of smooth solutions. For this numerical test, we consider a smooth initial data 

Po{x) = Pg + 0.2 sin(7rx). 



(6.J 



5' 



1 / TTln 

2 V Po 



which is initially at the local equilibrium zq = A (pQ,mQ). We apply periodic boundary condition in space 
and a CFL number A = 0.9. 

On the one hand, we compute as in the previous case an estimate of the numerical error for different 
meshes in space and time in order to evaluate the order of accuracy of the scheme with respect to h and 
for different regimes. In Figures [21 we present the plots in log scale of the error estimate given by 



f 



'I ~ P2h\\p + 



m 



2h\\p 



+ 



Z. 



2/illP' 



where h = {Ax, At) is the discretization parameter and || • ||p is the discrete norm with p G {1, 2, oo}. 
We observe that in this case the order of accuracy of our method is relatively closed to one both for 
and -L°° norms and the numerical error is not affected by the large variations of the relaxation parameter 

£. 
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Figure 2. Test 1 : Approximation of smooth solutions, (a) error and (b) L°° error 
at time t = 1 for different regimes with initial data 16. 8\) . 



On the other hand, we investigate numerically the long-time behavior of the solution to the Broadwell 
system ()6.4p . We consider / = (/-|_, /q, /-)^ a smooth solution to ()6.4p with periodic boundary conditions 
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in space. This solution is expected to converge to the global equilibrium {pg,mg) as t goes to +00. In 
order to observe this trend to equilibrium, we investigate the behavior of the quantities 

•Spit) = \\pit) - PgWi, Smit) = \\m{t) - nigWi, 

where the initial data is a perturbation of this global equilibrium is pg = 1, mg = and Zg = A {pg, rUg). 
In Figure El we present the time evolution of Sp{t), Sm{t) for the relaxation model and the hydrodynamic 
limit ()6.7p . Thus, we observe that the frequency of oscillations does not depend on e, whereas the 
damping rate is proportional to e. In other words, it appears that the trend to equilibrium is much 
faster (exponentially fast) in the rarefied regime, which corresponds to large values of e > than in the 
hydrodynamic one. In the hydrodynamic regime e = 0, the solution does not converge to the equilibrium 
since dissipative terms disappear. 




5 10 15 20 5 10 15 20 

Figure 3. Test 1 : Approximation of smooth solutions. Influence of the relaxation pa- 
rameter e: evolution of Sp (left) and Sm (Tight) for e = 0.5, 0.05, 0.005 and comparison 
with the hydrodynamic limit \6. 7| j (dashed black line). 



Test 2 : The Riemann problem. We now consider a discontinuous initial datum and propose numerical 
simulations to (j6.4p with different relaxation parameters e, from rarefied to fluid regimes. The initial 
condition is given by the local equilibrium 



(6.9) {po,mo,zo) 



(1,0,1), if-l<x<0, 
(0.25,0,0.125), ifO<3;<l. 



We approximate the Broadwell system in (—1, 1), with reflecting boundary condition and a CFL number 
A = 0.9 and present the evolution of the numerical solution. We only use 100 grid points, so that the time 
step is fixed equal to Ax = 0.002, and compare the results with the ones obtained from a fully explicit 
solver for which the time step has to be of the order of e for different values of the relaxation parameter. 

We first compute an approximation of the error with different meshes in space and time in order 
to evaluate the order of accuracy. In Figure HI we present the plots of £i{h) and £2ih) in log scale where 
h = (Ax, At) is the discretization parameter. In that case, the order of accuracy is not fixed and depends 
on the different regimes, but in the worst case the numerical error £i{h) is of order as in the estimate 
dSH]). 

In Figure [Sj we take e = 0.5 and 0.1. For such values of e, the problem is not stiff and this test is 
performed to compare the accuracy of our AP scheme with a fully explicit scheme (global Lax-Friedrichs 
method with slope limiters and explicit Euler discretization in time). The density p, the momentum 
(m, z), and the deviation to the equilibrium z — A{p, m) are plotted at different times t = 0.05, 0.2, 0.35 
and 0.5. At the kinetic regime, we can observe that our method gives the same accuracy as a standard 
fully explicit scheme. 
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Figure 4. Test 2 : The Riemann problem, (a) error and (h) error at time t = 1 
for different regimes, with initial data 16. 9\) . 
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Figure 5. Test 2 : The Riemann problem. Solution to the AP scheme (crosses) and 
solution the explicit solver (lines) with initial data \6.9(l in the kinetic regime: e = 0.5 for 
the density p, the momentum {m,z) and the deviation to the local equilibrium z — A{p,m). 



Next we investigate the cases of small values of e. The same time step for the AP scheme is used, 
whereas the fully explicit scheme requires it to be of order 0(e). We report the numerical results for 
e = 10~^ and e = 10~^ in Figure [H In this case, we add in the comparison the numerical solution to the 
hydrodynamic limit problem (16. 7p . obtained with a standard first order finite volume scheme. 

We observe that the AP scheme and the fully explicit scheme still agree, even if the time step is at 
least ten times larger with our method. Moreover, now that we are closer to the fluid regime, we see that 
the macroscopic quantities are in good agreement with the ones obtained with the limit problem. Yet 
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Figure 6. Test 2 : The Riemann problem. Solution to the AP scheme (crosses) and 
solution the explicit solver (lines) with initial data \6.9) in the kinetic regime: e = 0.1 for 
the density p, the momentum {m,z) and the deviation to the local equilibrium z — A{p,m). 



some differences between tlie AP and the explicit schemes can be observed for e = 10 this comes from 
the fact that we used a very small number of points for both discretizations. 

7. Conclusion 

In this paper we proposed a rigorous convergence proof of an asymptotic preserving numerical scheme 
applied to a system of transport equations with a nonlinear and stiff source term, for which the asymptotic 
limit is given by a conservation laws. We have proved the convergence of the approximate solution vfj 
to a nonlinear relaxation system, where e > is a physical parameter and h represents the discretization 
parameter. Uniform convergence with respect to e and h is proved and error estimates are also obtained. 
This theoretical result allows justify the efficiency of this new scheme for multi-scale problems. 
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